Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 17 August 2012 (MN KTeX style file v2.2) 



(N 



Of 



General relativistic radiation hydrodynamics of accretion flows: II. 
Treating stiff source terms and exploring physical limitations 

C. Roedig* 1 , O. Zanotti 2 , D. Alic 1 

1 Max-Planck-Institut fur Gravitationsphysik, Albert Einstein lnstitut, Am Muhlenberg 1, 14476 Golm, Germany 

2 University di Trento, Laboratorio di Matematica Applicata, Via Messiano 77, 1-38100 Trento, Italy 

17 August 2012 

ABSTRACT 

We present the implementation of an implicit-explicit (IMEX) Runge-Kutta numerical 
scheme for general relativistic hydrodynamics coupled to an optically thick radiation field 
in two existing GR-(magneto)hydrodynamics codes. We argue that the necessity of such an 
improvement arises naturally in most astrophysically relevant regimes where the optical thick- 
ness is high as the equations become stiff. By performing several simple one dimensional tests 
we verify the codes' new ability to deal with this stiffness and show consistency. Then, still in 
one spatial dimension, we compute a luminosity versus accretion rate diagram for the setup 
of spherical accretion onto a Schwarzschild black hole and find good agreement with previ- 
ous work which included more radiation processes than we currently have available. Lastly, 
we revisit the supersonic Bondi Hoyle Lyttleton (BHL) accretion in two dimensions where 
we can now present simulations of realistic temperatures, down to T ~ 10 6 K or less. Here 
we find that radiation pressure plays an important role, but also that these highly dynami- 
cal set-ups push our approximate treatment towards the limit of physical applicability. The 
main features of radiation hydrodynamics BHL flows manifest as (i) an effective adiabatic in- 
dex approaching 7 c ff ^4/3; (ii) accretion rates two orders of magnitude lower than without 
radiation pressure, but still super-Eddington; (iii) luminosity estimates around the Edding- 
ton limit, hence with an overall radiative efficiency as small as r] erLC ~ 10~ 2 ; (iv) strong 
departures from thermal equilibrium in shocked regions; (v) no appearance of the flip-flop in- 
stability. We conclude that the current optically thick approximation to the radiation transfer 
does give physically substantial improvements over the pure hydro also in set-ups departing 
from equilibrium, and, once accompanied by an optically thin treatment, is likely to provide a 
fundamental tool for investigating accretion flows in a large variety of astrophysical systems. 
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1 INTRODUCTION 

The field of numerical relativistic hydrodynamics has recently seen 
much progress in treating astrophysical systems under more and 
more realistic conditions. Because of the large computational costs 
involved, the inclusion of multi-dimensional general relativistic 
radiation hydrodynamics (GR-RHD) has been postponed for a long 
time, with the remarkable exception of neutrino transpo rt in the 
context of supernovae simulations fsee lLentz et alj J2012T) and ref- 
erences therein]. However, due to the increasing power of super- 
computers, the situation has started changing significantly in the 
last few years, and the inclusion of a photon-field is no longer re- 
garded as a remote possibility. 

This delay has, however, not been due to the fact that dynam- 
ical radiation fields are not regarded as a main ingredient, rather 
it is the inherent difficulty of solving the radiation transfer equa- 
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tiorQ. The cooling time-scales of a dynamical fluid may easily vary 
over several orders of magnitude within the computational domain. 
This then leads to characteristic propagation speeds for the pho- 
tons in optically thin regions that are much higher than the cou- 
pled fluid/photon speeds in optically thick regions. Not only are 
time-scales vastly different, but also additional spatial resolution 
is required whenever the coupling to the photon field induces small 
scale instabilities and turbulence. In addition, surfaces of astrophys- 
ical structures are typically not in local thermal equilibrium (LTE) 
and can cool very efficiently, usually on much shorter time scales 
than the dynamical ones. This problem becomes particularly se- 
vere when performing global simulations of astrophysical systems 
in which the principal force is gravity. In these cases, the spatial do- 



1 See lPomraninl jl973l) and lMihalas & Mihalaj Jl984l) fo r a com prehen- 
sive treatment of radiation hydrodynamics and lSchweized |l988j) for the 
extension to the relativistic case. 
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main must firstly be large enough to contain the entire astrophysical 
structure and secondly, it needs to resolve the influence of gravit)Q. 

Any such multi-scale problem is numerically extremely costly 
and it is thus important to formulate efficient algorithms that in- 
clude at least a leading order approximation to the various physics 
while still remaining computationally affordable. One of the most 
successful strategies was, and still is, represented by the so called 
projected s ymmet r ic trac e-free (PSTF) moment formalism intro- 
duced by iThornel d 198 lh . By defining moments of the radiation 
field similarly to how density, momentum and pressure of a fluid 
are defined as velocity moments of the corresponding distribution 
function, such a formalism provides an accurate, though still rea- 
sonably cheap, approximation to the solution of the radiation trans- 
fer equations. This approach is particularly appealing in the case 
of an optically thick medium, characterized by a strong coupling 
between matter and radiation. iFarris et af, (2008) were first to un- 
dertake the implementation of the corresponding radiation hydro- 
dynamics equations in a general relativisti c framework. A further 
step has been taken bv lShibata et alJi 201 lh. who ado pted the vari- 
able Eddington factor approach of Levermorel d 1984b to solve the 
relativistic radiation-hydrodynamics equations both in the optically 
thin and in the optically thick limit. This represents a significant 
progress with respect to simplified treatments, where effective cool- 
ing functions are introduced. 

In spite of all this progress, major numerical difficulties still 
prevent the application of such schemes to realistic astrophysical 
systems; one of them being t he pr esence of stiff source terms. 
For example, in lZanotti et al.l d201ll) (hereafter paper l), after im- 
pleme nting and testing the framework suggested by IFarris et all 
( 2008), we studied the Bondi Hoyle Lyttleton (BHL) accretion flow 
onto a black hole, but we could only treat unrealistically high fluid 
temperatures of the order of ~ 10 9 K or above. Though simpli- 
fied, the BHL flow can effectively help our understanding of those 
compact sources accreting matter with a reduced amount of angu- 
lar momentum, and is currently applied to the s tudy of both High 
Mass X-ray Binaries (Hadrava & Cechurall2012|) and of the merg- 
ing of supermassive black hole binaries fsee IPfeiffeJ 120121) and 
references therein]. 

In this paper, we address the problem of treating the optically 
thick regime compatible with the conservative formulation used 
in Eulerian GR-MHD codes, while at the same time coping with 
the stiffness of the source terms. A s a stiff solver, we choo se the 
implicit-explicit (IMEX) scheme by Pares chT& Russol J20051) . im- 
plement it in both WHISK^Q and EchcQ, and test the codes against 
each other. As the two codes contain internal differences, such as 
scheduling and general infrastructure, it is very useful to validate 
them both at this stage, even though the main part of the simula- 
tions shown in this paper are performed with ECHO, because of its 
spherical, non-uniform gric0. 

The paper is organized as follows. In Sec. [2] we describe the 
treatment of the radiation stiff source terms. We detail an IMEX 
Runge Kutta scheme as our time integration stiff-solver. Sec. [3] 
presents the verification of our new scheme through a selected sam- 
ple of stiff shock tube problems. Turning towards astrophysical ap- 
plications, we first present in Sec.|4]the results for spherical accre- 



A complementary approach, which is not covered here, is to model not 
a global system, but only a small, representative region e.g. a shearing box. 
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5 WHISKY uses Cartesian adaptive mesh refinement, which is less suited 
for spherical models. 



tion in a regime that was constructed to be particularly challenging 
for the numerics. We also present a physical Michel solution and 
compare it with previous results. Abandoning spherical symmetry, 
we devote Sec. [5] to the study of the radiation hydrodynamics of 
BHL accretion in two dimensions. Finally, in Sec.|6]we offer a brief 
summary and our conclusions. 

Throughout the paper, we set the speed of light c = 1, and the 
gravitational constant G to a pure number. We extend the geometric 
units by setting m p /kB = 1, where m p is the mass of the proton, 
while ks is the Boltzmann constant. However, we have maintained 
c, G, and ks in a explicit form in those expressions of particular 
physical interest. We refer the interested reader to Appendix lAl of 
paperl for the system of extended geometrized units. 



2 RADIATION HYDRODYNAMICS IN THE STIFF 
REGIME 

2.1 Formulation of the GR-RHD equations 

In this section, we first review the set of equations that we use to ap- 
proximate general relativis tic radiation -hydrodynamics in the dif- 
fusion limit, as derived in Farri s et al] d2008l) and already imple- 
mented and verified in paperl. The properties of the fluid immersed 
in the radiation field are described by the momentum-energy tensor, 
which is given by 



and comprises a matter contribution 



a/3 



= phu a u p + Pg 



a/3 



and a radiation contribution 

af3 _ 1 



T, 



I„N a N p dvdn, 



(1) 



(2) 



(3) 



where g afi is the metric of the spacetime, u a is the four-velocity of 
the fluid, p, h = 1 + e + P/p, e, and P are the rest-mass density, 
the specific enthalpy, the specific internal energy, and the thermal 
pressure, respectively, while I v — I v (x a , N r , v) is the specific in- 
tensity of the radiation. We note that N a defines propagation direc- 
tion of the photon with frequency v, while dQ, is the infinitesimal 
solid angle around N a . All of these quantities are measured in the 
comoving frame of the fluid. The thermal pressure is related to p 
and e through an equation of state (EoS), which we take to be that 
of the ideal-gas, with constant adiabatic index 7, i.e. 



P = pe(j - 1) . 



(4) 



In terms of the moments of the radiation field jThornd fl98lh . 
the radiation energy-m omentum tensor T"^ can be rewritten as 
fasten & Spiegelll 19761) 



Tf = (E T + Vr)u a u p + F r a u p + u a F p + Vrff 



(5) 



where E r and P r are the radiation energy density and pressure, 
respectively. We make the additional assumption that the radiation 
field is approximately isotropic, in the sense that Vi = E T /3, while 
the radiation flux is not constrained to zero, but is allowed to take 
small values such that / E T <C 1. Thus the equations governing 
the evolution of the system are: 



V a {pu a ) =0, 
V Q T Q ' 9 = 0, 
V a T r afi 



-G 



(6) 
(7) 
(8) 
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where G" = G™ (I, x* , X s )' called the radiation four-force density, 
depends on the specific intensity and on the opacities of the matter 
interaction. As in paperl, we drop all frequency dependencies and 
allow for small deviations from LTE. We consider bremsstrahlung 
and Thomson scattering (i.e. x an d x") as processes of absorption 
and scattering. Using the Planck function, B, it is t hen possible 
to wr ite the radiation four-force in covariant form as dFarris et"al] 
2008) 



paperl for more details): 



(9) 



In Eq. (|9) we have introduced the equilibrium black-body intensity 
AtvB = aradTfluid, where Tfl u id is the temperature of the fluid and 
a ra d is the radiation constant. We estimate the temperature from the 
ideal-gas EoS via the expression 



Tfluid 



m p P 
k B p 



(10) 



where, ks is the Boltzmann constant and m v the rest-mass of the 
proton. We stress that the method allows for deviations from ther- 
mal equilibrium, namely with E r 7^ 4irB. As shown in paperl , 
after adopting the 3 + 1 split of spacetime dArnowitt et al 1 ll962h 
the GR-RHD equations can be written in conservative form as 



d t U + d,T l = S , 



(11) 



where the vector of conserved variables U and the fluxes J 71 are 
given by 



D 



Sj 

W = \ ~ u 

u t 

while the sources are 



av'D- 


- p i D 


aW)- 


- P% 


aS l - 


/3'U 


aSl — 




a(7? r )J - 





(12) 



\aW lk dj^ k + Sidjp* - UdjOt + a{G T )j 

\W ik ftd^ ih + Wjd j p i - S ] dja + a 2 Gl 

iR^P'd^k + {Rtfidjp - Stdja ~ a 2 Gt 

\aRi k djMk + (Sthdjtf - Urdja - a(G r ) d 

(13) 

We note that a, j3, and k are the lapse, the shift, and the deter- 
minant of the spatial metric, respectively, while v x and Y are the 
three-velocity and the Lorentz factor of the fluid with respect to 
the Eulerian observer. In the Eqs. l !12t and J X 3 b several more terms 
have been defined, which we report below for completeness (c.f. 





= phY 2 v l v ] + Pi' j , 


(14) 


S' 


= phY 2 v l , 


(15) 


u 


= phY 2 - P , 


(16) 




= -jsyrW + Y(fiv j + fitf) + v t iP , 


(17) 


si 


= i J e r rV + r(aF r V + /;), 


(18) 


Ur 


= ^E T Y 2 + 2aYFi - ^ , 


(19) 


Fi 


ViF* Vif* 
a — f3iV l a 


(20) 



2.2 



Description of the IMEX scheme for radiation 
hydrodynamics 



2.2.1 General concepts 

A relevant feature of the radiation hydrodynamics equations 1 11 It 
is that they contain sources for the radiation field that may easily 
become stiff, depending on the physical conditions under consid- 
eration. When stiffness is treated by resorting to implicit-explicit 
(IMEX) Runge-Kutta (RK) scheme^ it is important to split the 
conservative variables U in two subsets {X, Y}, with {-X"} con- 
taining the variables that are affected by stiffness, and {Y} con- 
taining those that are not. IMEX Runge-Kutta methods are based 
on an implicit discretisation for the stiff terms and on an explicit 
one for the non-stiff ter ms. They have been exten sively discussed 
in a series of papers by IPareschi & Russol J2005t) . and some re- 
cent applications have been presented in special relativistic re- 
sistive MHD by Palenzuela et al. (2 0"5^), in general relativistic 
force-free electrodynamics bvlAlic et al.l i2012|) and in genera l rel- 
ati vistic resistive MHD by Bucciantini & Del Zannal j2012h and 



by iDio nvsopoul ou et al. I J2012h . In full generality, the hyperbolic 
equations for the two sets of variables {X , Y} are split as 



d t Y 
d t X 



F Y (X,Y) 
Fx(X,Y) 



Rx(X,Y), 



(21) 
(22) 



where the operator Fy contains both the first spatial derivatives of 
Y and non-stiff source terms, the operator Fx contains both the 
first spatial derivatives of X and non-stiff source terms, while the 
operator Rx contains the stiff source terms affecting the variables 
X. Each Runge-Kutta sub-stage of the IMEX scheme can be di- 
vided in two parts. 

(i) In the first part, the explicit intermediate values 
{X*' 1 , Y*' 1 } of each sub-stage i are computed as 

i-l 

Y*'' = Y n + AtY / a l3 F Y [U u) ], (23) 

3=1 

i-l i-l 

X*' 1 = X n + AtY,a tJ F x [U u) ]+AtY,a tj Rx[U u) ], 
3=1 3=1 

(24) 

where one might note that the summation stops at (i — 1), in or- 
der to avoid the appearance of the implicit terms at this stage. The 

6 An alternative approach to solve the special relativistic RHD equations 
in a moderately stiff reg ime has been consid ered in one-dimensional La- 
grangian simulations by Dumbser et al] 120121) . 
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matrices and (ay) are v x v square matrices. In this paper, 
we use v = 4 (see also Appendix [B}, whereas, in general, the ma- 
trix coe fficients and dimensions change with the desired number of 
staged jPareschi & Russoll2005t) . 

(ii) In the second part, the non-stiff variables are directly ad- 
vanced to the status of sub-stage Runge-Kutta variables, namely 

y W = Y *>* , (25) 

while the stiff variables need to be corrected as 

x d) = m{Y*- 1 ) [x*'* + auAtKx(Y**)] . (26) 

The vector Kx (Y) on the right hand side of Eq. J26l >, which does 
not depend on the stiff variables X, results from the decomposition 
of R X (X,Y) as 

Rx(X,Y) = A{Y)X + K X {Y), (27) 

while the matrix M is given by jPalenzuelaetal] l2009fl 

M(Y*'') = [/ - a rr AtA(Y*^ ~* , (28) 
where I is the identity matrix. 

For each RK sub-stage, {X^\ Y^} is computed as described 
above, and finally the time-update is performed as 

V V 

U n+1 =U n + At^2mF[U^] + M^2wiR[U (i> ] , (29) 

where Wi and Wj are coefficient vectors. In most of the applications 
presented in this paper, we have adopted the SSP3(4, 3, 3) (Strong 
Stability Preserving of order three) IMEX Runge-Kutta scheme. 
The notation SSPk(s,<r, p) is adopted to specify the order of the 
SSP scheme (fc), the number of stages of the implicit scheme (s), 
the number of stages of the explicit scheme (a), and the order of the 



IMEX scheme (p) jPareschi & Ru sso 2005). The coefficient tables 
employed in this paper are listed in the Appendix |B1 

2.2.2 Specification to radiation hydrodynamics 

Because of the complexity of the GR-RHD equations, isolating the 
term (or the terms) that are responsible for the stiffness is not a triv- 
ial task, although we can certainly say that such terms are contained 
in the radiation four-force G" . According to the logic of the IMEX 
scheme just described, we identify {X} with the radiation hydro- 
dynamical variables {U T , {S T )j} that are affected by stiffness, and 
{Y} with {D, Sj , U}, that remain unaffected. 

As highlighted above, the IMEX scheme requires the stiff 
source terms Rx to be decomposed according to Eq. j27\ . We 
therefore write the radiation four-force G" in terms of the con- 
servative variables of the radiation field. To this extent, we rewrite 
Eq. d 1 8b and Eq. dl9t to find the radiation energy density E T and 
the fluxes F" in terms of U r and (S r )i as 

E T = -3T 2 W [2(S r ) k v k + (7 r (l/r 2 - 2)] , (30) 
Fi = ^-W^4U r {T 2 -1) + {4F 2 -l)(Sr) k v k ] , (31) 
(/*)< = ^-^Tvi-aiFrfvi, (32) 

7 Note that the global order of an IMEX scheme does not uniquely deter- 
mine the number of sub-stages. 

8 We stress that the form of M given by Eq. d28t is only valid for the 
decomposition as done in Eq. (27). 

where W = 1/(1 + 2r 2 ). In this way, and after some simple alge- 
bra, we can rewrite the radiation four force as 



G\ = - | [x'a r Tj 4 nid + U r (2 X s (l - 3W) - X ) + (5 r )fc« fc (x + x' W - 2))] , (33) 

(Gr)i = - xVlti^r + l^±*12(S r)l + UrTvi [ X *(l - 4W) + 2 X S (W - 1)] + (S^Fv, [ X t (2W - 1) + X s (2 - W)] . 

(34) 



We note that the right hand sides of d33t and d34t do not con- 
tain the set of variables {Y}, while they do contain the conserved 
variables {X}, which always appear with a multiplication factor 
containing either x or x" ■ This is an indication that, depending 
on the values assumed by the opacities, such source terms may be- 
come stiff, but only for the radiation variables. This means that 
the vector of sources given by Eq. J 1 3t will be split in two parts, 



S = <S e + Si. The first one, 



\aW ih d^ lk + Sidjp* - Ud 3 a + a(G t ) s 

S e = V^ | W ik P j djr ik + Wi j d 3 P l -S J d 3 a + a 2 Gl 

iR^ftdjMk + {Rr) i i d j p i - Sid.a 

laRfdjMk + (Sr^djP - U r d 3 a 

(35) 

will be absorbed into the operators Fy and Fx m Eqs. d2U 
and d22t . because it does not contain stiff terms. The second part, 
on the other hand, which contains the genuinely stiff terms for the 
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radiation variables {X}, is 



Si = yft 








(36) 



-a(G t )j 

and its non-zero components are identified with i?x(-X^,^) in 
Eq. {22j. After using Eq. i33\ and Eq. l !34t . it is possible to further 
decompose Rx as prescribed by Eq. j27\ as 



—a G r 
















+ 




-a(G t )i . 











(37) 

where the coefficients of the matrix A(y) are specified in the 
Appendix [B] The components of the vector K x (the second term 
on the right hand side of Eq. ( 137b ) do not depend on the stiff 
variables X, but only on the temperature Tfl u id. We note that, in 
the actual implementation of the IMEX Runge-Kutta scheme, the 
correction to the implicit variables X' 1 ' dictated by Eq. d26t is 
performed when the conversion from the conservative variables U 
to the primitive variables is performed. 



2.3 Numerical tools 

For reasons of flexibility, cross-verification and in view of future 
projects, we have implemented the GR-RHD equations in their 
IMEX version in two different numerical codes. 

The first one is a modification of the WHISKY code, which 
implements the general relativisti c resistive magnetohydrod ynam- 
ics formalism WhiskyRMHD dDionvsopoulou et alj|2012h . We 
use the numerical methods provided by the original WHISKY 
code documented in (Baio ttiet alj|2003l ; iGiacomazzo & Rezzollal 
120071) . namely an HLLE approximate Riemann solver and a sec- 
ond order TVD slope limiter method for the reconstruction of 
the primitives. The infrastructure as well as the solution of the 
Einstein equations is prov ided by the Cactus Computational 
Toolkit dLoffler et al1l2012h . The implementation of the GR-RHD 
equations in WhiskyRMHD required modifications mainly in the 
sources and the routine which recovers the primitives from the con- 
servative variables. In order to deal with the stiffness of the source 
terms, we have modified the MoL thorn (part of Einstein ToolkiQ), 
by including second and third-order IMEX Runge-Kutta time inte- 
grators. 



methods. Time integration is possible in either second or third- 
order IMEX Runge-Kutta. Previously in paperl, ECHO had been 
extended to allow the solution of the non-stiff GR-RHD equations 
in the optically thick regime. 

In both WHISKY and ECHO, our implementation of the stiff 
GR-RHD equations does not allow for a treatment of the optically 
thin regime. Therefore, all the tests and applications described in 
this paper are limited to the optically thick regime, while we post- 
pone an accurate analysis of the variable Eddington factor approach 
to a future work. 

Finally, we note that the increase of computational cost when 
changing from an explicit RK of order k to a RK-IMEX of the same 
order k is approximately given by the ratio of the number of sub- 
stages required by the IMEX to the number of sub-stages required 
by the explicit RK. For the SSP3-IMEX scheme, compared to the 
explicit RK3, such a nominal ratio is given by 5/3 ~ 1.67 and 
in both our implementations we have measured an effective factor 
~ 1.8 increase. 



3 VERIFICATION OF THE SCHEME 

In paperl we had considered a number of shock-tube tests in 
which nonlinear radiation-hydrodynamic waves propagate. The 
semi-analytic solution that is used for comparison wit h the numer- 
ical on e has been obtained following the strategy of iFarris et al.l 
(2008), and it requires the solution of the following system of ordi- 
nary differential equations 



d x V(P) = S(P) : 



(38) 



where 



( t \ 



P 
P 

u* 
E, 
F, 



U = 



pu 

rj-iOx 
rpXX 



t: 



•Ox 



s = 







-G r ° 
V -Gr ) 



Ui, lli and Uz are constant in x, while only T r and T T XX need to 
be solved for. These tests can be used to monitor the ability of the 
code to deal with the stiff regime, by simply increasing the thermal 
opacity «* (the scattering opacity n s g is set to zero). When this is 
done, the semi-analytic solution of the ODE system ([38} can be ob- 
tained with an ODE solver for stiff systems iPress et alj 19921) . The 
initial states of the two tests that we have considered are reported in 
Table[T]and are chosen in such a way that the discontinuity front at 
x = remains stationary, namely it is comoving with the Eulerian 
observer. LTE is assumed at both ends x = ±X, with X — 20, and 
this is obtained by adopting a fictitious value of the radiation con- 
stant a rac i, namely a ra d = E t ,l/TI, which is then used to compute 
Er.R = a ra dTj| (here the indices L and R indicate the "left" and 
"right" states, respectively). We note that tests No. 1 and 2 in Ta- 
ble [p a re the same of tests No. 3 and 4 in Table [T1 of IZanotti et al.l 
( 120 1 lh . apart for the value of «*, which controls the stiffness of the 
problem. After setting 800 grid-points in the x direction, we have 
increased the value of t^ g to the maximum value affordable by the 
numerical scheme, while keeping the Ccfl parameter unchanged 
and equal to 0.25. For example, K g has been increased from 0.3 

9 jThe Einstein Toolkit: A Community Computational Infrastructure for Relatjvist^^ fMVM"- Na 1 - and from °- 08 to °' 7 in test No - 2 - Each test 
i cac | ^ is evolved in time until stationarity is reached, and the results are 

10 See also Bucciantini & Del Zanm] 4201 lh for a recent extension of 
ECHO to dynamical spacetimes within the conformally flat approximation. 



The second code is based on ECHO dDel Zanna et~aT1l2007l) . 
which provides a numerical platform for the solution of the 
GRMHD equations in stationary background spacetimef°l. It em- 
ploys a high-order shock-capturing Godunov scheme with a two- 
waves HLL Riemann solver, while the spatial reconstruction of 
the primitive variables can be obtained by linear and non-linear 



shown in Figure[T] where the numerical solution is compared to the 
semi-analytic one in the two cases considered. 
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Table 1. DESCRIPTION OF THE INITIAL DATA - in the shock-tube tests with radiation field. The different columns refer respectively to: the test considered, 
the adiabatic index, the radiation constant and the thermal opacity. Also reported are the rest-mass density, pressure, velocity and radiation energy density in 
the "left" (L) and "right" (R) states. 
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Figure 1. SHOCK TUBES - Solution of the test No. 1 (left panel) and No. 2 (right panel). From top to bottom the panels report the rest-mass density, the 
velocity and the radiation energy density. In both cases 800 grid-points have been used with Ccfl = 0.25 and RKIMEX2. The tests are performed with the 
WHISKY code, employing TVD reconstruction and minmod limiter. 



It should also be noted that shock tube problems do not rep- 
resent an ideal set-up to highlight the ability of the scheme in han- 
dling the stiffness of the source terms, since strong discontinuities 
are by themselves a challenge for any numerical method. As a re- 
sult, we have performed an additional and peculiar shock tube prob- 
lem, test No. 3, which has equal left and right states, except for the 
velocity. In this case, two shock waves propagate in opposite di- 
rection, no stationary solution is obtained, but a much higher value 
of K g can be used, namely k* = 1000. Figure [2] reports the cor- 
responding solution at time t — 15, and also shows the very good 
agreement between the results obtained with WHISKY and ECHO. 



4 SPHERICAL ACCRETION 

Having introduced the numerical tools for the treatment of the stiff 
source terms typical of GR-RHD, we now focus on a problem 
that has been the subject of several astrophysical analyses, namely 
spherical accretion onto a black hole. In the first part of this §, we 
present an additional test of our numerical scheme, brought in the 
stiff regime by assuming unphysically large cross-sections. On the 
other hand, in the sec ond part, we c hoose physical parameters to 
model the solution by Michel] J 1972t) in an astrophysical context. 

Transonic accretion onto a non-rotating black hole in the 
presence of an isotropic radiation field has been studied in great 
detail by several research groups over the years. In the opti- 
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Figure 2. SHOCK TUBES - Solution of the test No. 3 at time t = 15. 
From top to bottom the panels report the rest-mass density, the velocity and 
the radiation energy density. In both cases 800 grid-points have been used 
with Cqfl = 0.25 and RKIMEX2. The tests are performed with both the 
WHISKY and ECHO codes, employing TVD reconstruction and MC limiter. 
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cally thick regime, the stationary solution was investigated under 
different approximations and by foc u sing on different e mission 
mechanisms b y|Maraschi et al. ( 1974), Kafka & Meszaros ( 1976), 
Vitelld dl978|),lGillman & Stellingwerj dl98CMFlammang| dl982l) . 
Nobili et all dl99l | ). The tim e dep endent solution , was c onsidered 
by lGilden& Wheeled dl980h and IZampieri et al.1 dl996t) . The lat- 
ter, in particular, solved via a Lagrangian code the radiation trans- 
fer equations using the PSTF moment formalism^ truncated at 
the first two moment equations. Because of the limiting approx- 
imations assumed, and in particular because of the lack of Comp- 
tonization effects, our analysis should not be regarded as an attempt 
to improve with respect to the above mentioned works, but rather as 
a preliminary study in view of further developments. We also note 
that multidimensional simulations with an Eulerian code have been 
recently performed by Fra gile et al.l d2012l) obtaining promising re- 
sults. 

Our initial conditions a re given by th e fluid spherically sym- 
metric transonic solution of Michel (1972), which is stationary in 
the absence of a radiation field. The free parameters of the fluid 
solution are the critical radius r c and the rest mass density at the 
critical radius p c . We choose a black hole with mass M — 2.5Mq , 
while the adiabatic index of the fluid is 7 = 4/3. The initial radi- 
ation field is initialized to a negligible energy density, while radia- 
tion fluxes are set to zero. As a first test, aimed at showing the abil- 
ity of the numerical scheme in handling the stiff regime, we have 
considered an unphysical setup with p c = 0.02, r c = 8.0, and a 
high uniform value of the thermal opacity, K g = 10 lj . The test is 
performed in Boyer-Lindquist coordinates with 2.5 < r < 200 us- 
ing N — 300 radial grid points. The SSP3-IMEX scheme has been 
adopted, with the MC limiter for the spatial reconstruction. Figure[3] 
shows the profiles of the rest-mass density, the radial velocity and 
the radiation energy density (from top to bottom) at time t = 1000. 
We stress that, if the IMEX scheme is not available, and the evo- 
lution is performed through a fully explicit Runge-Kutta scheme, 
this test can be successfully repeated at the same Ccfl = 0.2 only 
with a value of Kg<l. 

Having done that, we have concentrated on a sequence of more 
realistic models, all of them with p c = 9.88 x 10~ 9 cgs, but with 
different critical radii, chosen in the range between r c — 800 and 
r c — 7000, in order to control the accretion rate. The test is per- 
formed in Kerr-Schild coordinates with 1.0 < r < 1000 using 
N — 3200 radial grid points The evolution is stopped when sta- 
tionarity in the L2~ norms of all of the variables has been reached, 
which may require a final time as long as t — 400000 in code units. 
The scheme employed is the SSP3-IMEX, with MC limiter. 

Special attention has to be paid to the boundary conditions 
at the outer radial grid poin t, for which we hav e followed closely 
the discussion presented bv lNobili et all dl99ll) . In particular, ze- 
roth order extrapolation (copy of variables) is adopted for the gas 
pressure and for the density. This guarantees that the tempera- 
ture has zero gradient. At the same time we want to make sure 
that at large radii the radiation field streams radially, namely that 
E oc f r oc r~ 2 . This translates into the condition 

dinr 

which can be easily implemented. Finally, we fix the accretion rate 
at the outer boundary to the value possessed by the initial configu- 
ration. At the inner radial boundary, on the other hand, zeroth order 

11 The first time dependent problem s adopting the PSTF formalism were 
presented in Rezzolla & Milled Jl994[) . 




Figure 3. STIFF SPHERICAL ACCRETION - Numerical solution at time 
t = 1000. From top to bottom the panels report the rest-mass density, 
the velocity, the radiation energy density. An artificial k* = 1.0 X 10 
has been adopted to highlight the ability of the code to treat the stiff regime. 
N r = 300 grid-points have been used with Ccfl = 0.2, MC reconstruc- 
tion and SSP3-IMEX. 



extrapolation is adopted for all of the variables. The evolution is 
performed considering both the contribution of the bremsstrahlung 
opacity and of the Thomson scattering opacity for electrons. We 
note that during the evolution the radiation flux remains typically 
two orders of magnitude smaller than the radiation energy, thus 
maintaining the code in the physical regime for which it was de- 
signed. After an initial relaxation, the system converges to a differ- 
ent stationary configuration characterized by a non-zero radiation 
flux. The solution is optically thick in all the models for r ^ 100, 
while it becomes marginally optically thick at large radii. From the 
radiation flux we compute the luminosity as L = Aivr 2 f r . 

Fig. [4] reports the results of our simulation tests in the dia- 
gram (M/MEdd, L/IjEdd), where the luminosity is computed at 
r — 200. Although o ur data resemble t he high luminosity branch 
reported in Fig. 1 bv lNobili et all dl 99 ll) . a close comparison with 
their results is not really possible, since Comptonization, bound- 
bound transitions and free-bound transitions are not taken into ac- 
count in our analysis. In particular, the absence of pre-heating ef- 
fects does not allow us to verify the onset of strong thermal insta- 
bilities producing hyd rodynamic sho c k wav es that propagate out- 
ward, as reported by IZampieri et al.l d 19961) . In spite of this, the 
test we have performed is very relevant. In fact, by using a n en- 
tirely different proced ure with respect to iNobili et al.l dl99 ll) and 
IZampieri et at] (1996), it confirms the existence of a high lumi- 
nosity branch in the diagram (M /A^Edd, L/L^dd), which corre- 
sponds to the optically thick regime. A more extended analysis of 
this problem, by including additional contributions to the opacity, a 
treatment of the Comptonization and the effect of a spinning black 
hole will be the focus of a separate and dedicated work. 



5 BONDI-HOYLE-LYTTLETON (BHL) ACCRETION 

This Section deals with the application of our new scheme to simple 
astrophysical models departing from spherical symmetry. We re- 



8 C. Roedig, O. Zanotti, D. Alic 



io- 2 



j 10- 3 r 



10" 




/.if Nobili el al. (91) 



10° 



10 1 

M/M Edd 



10 2 




2xl0 4 



Figure 4. Spherical accretion: luminosity(accretion rate M) 
- in Eddington units. Luminosity L was extracted either at constant opti- 
cal depth t or at constant radius r. Additiona lly in red, we show the high 
luminosity branch found in iNobiliet alj|l99ll) for comparison. 



visit the BHL accretion flow, that we already described in some de- 
tail in paperl, and whose initial conditions are briefly summarized 
in Sec. 15. II After showing consistency with the non-stiff solver, we 
illustrate the effectiveness of the IMEX by treating models of low 
temperature, which is the key parameter responsible for numerical 
difficulty. Only now, it becomes feasible to treat astrophysical tem- 
peratures that are few orders of magnitude lower than in paperl and 
as astrophysically realistic as our approach can allow at this stage 
(c.f. conclusion for more discussion). Having thus reached the lim- 
its imposed by the physical assumptions of the current treatment, 
we now analyse the dynamics of the fluid, the occurrence of shocks 
and the possible observational quantities, with particular attention 
to the computation of the luminosity. We stress that it is not our in- 
tention to perform a systematic analysis of the full parameter space. 



Figure 5. Consistency test of the new IMEX scheme - Time evolu- 
tion of a perturbed BHL model with v ia { = 0.18 and c s>00 = 0.07 (model 
p.V18.cs07 of paperl) with the IMEX (solid line) and with the non-IMEX 
version of the code (dashed line). The lower panel shows the light curve, 
whereas the accretion rate is plotted in the top panel. All curves are shown 
in Eddington units. 



the two parameters Voc ± and c SjOO0 J3and a prefix denoting per- 
turbation (if applicable) in our naming scheme as V000.1 -Cs, 000.1 • 

The adiabatic index of the fluid is 7 = 5/3. 

The computational grid consists of N r x N$ numerical cells in 
the radial and angular directions, respectively, covering a compu- 
tational domain extending from r m m = 2.1 M to r max = 200 M 
and from m i n = to </> ma x = 2n. For our fiducial simulation we 
have chosen N r = 1536 and = 300, but have also verified that 
the results are not sensitive to the resolution used or to the location 
of the outer boundary. 



5.1 Initial conditions for the BHL accretion flow 

We perform two dimensional numerical simulations of a BH L ac- 
cretion flow faovle & Lvttletonlll939l ; lBondi & Hovlelll944l) onto 
a Schwarzschild black hole of galactic size with Mbh = 3.6 x 
1O 6 M0. The initial conditions considered are similar to those 
adopted in paperl, with a velo city field that is sp ecified in terms 
of an asymptotic velocity Uoo, jFont & Ibanez|[l998h 

v 1 — V if""«cx) cos <j) , (40) 
v <f> = _ vr^Woo sin0 , (41) 

where x 1 -' are the components of the 3-metric and (j> is the az- 
imuthal angle in Boyer-Lindquist coordinates. The radiation field 
is initialized to a uniform and small energy density E r , such that 
the radiation temperature r rad = (-Er/cirad) 1 ^ 4 w 1.5 x 10 5 i\". 
Additional free parameters are the asymptotic sound speed c aj00 , 
and the asymptotic pressure, from which the asymptotic rest-mass 
density follows directly. The resulting configuration relaxes to 
a different and stationary one, on a timescale that depends on the 
parameters chosen. Keeping to nomenclature of paperl, we encode 



5.2 Consistency test 

Before going to new models, we first carried out a consistency test 
using a representative model with Mach number Moo = 2.57 
(model p.V18.cs07 of paperl) and reproducing it with the present 
new IMEX- version of ECHO. As shown in Fig. [5] the IMEX ver- 
sion reproduces the light curves and the accretion rates obtained 
with the purely explicit version of the code. Moreover, by using the 
IMEX scheme, it is now possible to extend the evolution to later 
times, whereas the previous version of the code required reducing 
the Ccfl to values smaller than 0.01, making such long evolu- 
tions practically unfeasible. This test confirms that the new scheme 
is verified also in a non-trivial two-dimensional application and that 
the use of the IMEX offers clear advantages in terms of computa- 
tional resources. 



1 Here, subscripts 0.1 denote the normalisation in units of 0.1c, so 
"oo n ! = Utx>/(0.1 c). Therefore, the model V07.cs03 has Voa = 0.07 
and c s>00 = 0.03. 
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5.3 Results 

In the following we examine the behaviour of three models, with 
two different initial soundspeeds c s ,oo and the same asymptotic 
velocity v x . The prefix sp is used to denote "strongly perturbed" 
which means that the initial asymptotic pressure is lowered by two 
orders of magnitude with respect to the equilibrium value. This is 
done with the main purpose of producing models with even lower 
temperatures. 



Measured physical quantities In addition to the primitive vari- 
ables provided by the code, we calculate several physical quanti- 
ties: the accretion rate in Eddington units, M, the luminosity in 
Eddington units, L, the radiation equivalent temperature, T ra( j = 
(-Er/cirad) 1 , the fluid temperature, Tfl u id, the effective adiabatic 
index 7 cff : 



_ 5/2 + 20g + Wq 2 
7cff ~ (3/2 + 12g)(l + 9 ) 

and the local Mach number A4 : 

r 7 



with q = Vt/P ■ 



M 



with 




(42) 



(43) 



Moreover, as discussed in paperl, we define an effective BHL lumi- 
nosity efficiency t] erlL , that takes into account the injected energy 
at infinity as 

L 

Vbhc = — 7T-- — • C 44 ) 



Ma. 



= 2 + iMootii 



We measure several quantities Q as volume weighted averages 
over all grid elements i, thus defining the pointy brackets as 



(Q) = 



J2i ridrn 



(45) 



The rate of entropy generation is measured according to Eq. iA5i , 
which is an appropriate approximation for a coupled photon fluid 
plasma in a quasi-stationarity state. When we extract the luminos- 
ity, we choose a surface of constant optical depth r>10. We argue 
this is reasonable because only if r 3> 1 the system is still in a 
regime where the approximation with the diffusion limit is valid. 
This is also realistic, when thinkin g of actual observations, where 
measurements are taken at constant 13 ! r. For a discussion of how 
the luminosity estimates change with respect to paperl, see Ap- 
pendix [A] Suffice it to say that this luminosity is a tracer of the 
outwards radial fluxes and that different possible extractions agree 
within the current error bars. 

The optical depth r is computed in post processing as in paperl 



(X +X S )ds, 



(46) 



where we have assumed a constant characteristic lengthscale 
L — 10. All other quantities are standard and reported in cgs units. 
Selected results are shown in Figs.[6]to[8]and summarised in Tab.[2] 



In the case of stars, for instance, this is usually taken as r ~ 2/3. 
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Figure 7. Comparison of fluid and radiation temperatures - 
Volume weighted averages according to Eq. {45) for all three BHL models 
sp.V07.cs05, sp.V07.cs03 and V07.cs03 as a function of time. Tempera- 
tures are given in Kelvin. 



BHL dynamics dominated by radiation quantities First of all, 
we note that the qualitative dynamics of all BHL models is the same 
as described i n pape r l and can be summari s ed as follows [see also 
iPetrich et all dl989t) : iFont & Ibanezl dl998h: iDonmez et al.l J201lh 
for the hydrodynamics case and lPennen ( 120111) for the magnetohy- 
drodynamics one]. Initially, a narrow, hot shock cone forms down- 
stream of the accretor and the plasma is fluid-pressure dominated. 
Progressively, the radiation field builds up strength until the radia- 
tion pressure becomes similar to the fluid pressure. At this point, the 
shock cone becomes unstable, oscillating from one side of the ac- 
cretor to the other, until it finally reverses into the upstream domain 
as a bow shock. From now on the radiation pressure exceeds the 
fluid pressure, the effective adiabatic index approaches the value 
~ 4/3 and, at the same time, the density (and correspondingly the 
optical depth) decreases in most parts of the numerical domain. Af- 
ter the upstream shock has moved out of the numerical domain and 
expelled a significant amount of mass, a new, low density equilib- 
rium is formed in which there is a smaller shock cone in the down- 
stream region (this is illustrated well by Fig. 3 of paperl, which 
shows a comparison of BHL flows with and without the radiation 
field). 

The central improvement over paperl is shown in Fig. [6] where the 
two-dimensional maps of the optical depth (left) and of the fluid 
temperature (right) are shown for two different models, both of 
them at the final quasi-stationary state. The top panels, in partic- 
ular, show that for the strongly perturbed model sp.V07.cs03 large 
parts of the upstream region settle down to temperatures of the or- 
der Tfluid <10 6 K, a value which could not be reached before due 
to the stiffness of the equations. The corresponding unperturbed 
model, V07.cs03, shown in the bottom panels, has upstream tem- 
peratures as high as Tfl u id<5 x 10 9 K, while both of the models 
have significantly high optical thickness. It is important to note that 
at this stage of the evolution, namely after the "reversal" of the 
shock cone, the effective adiabatic index of all three models is very 
close to 7 c ff ~ 4/3, and therefore behaving like an effective pho- 
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Figure 6. 2D Optical depth and fluid temperature of perturbed model sp.V07.cs03 and unperturbed models V07.cs03- 
Both models are shown at stationary state: (Top) t = 7.71 X 10 4 M for model sp.V07.cs03 and (Bottom) t = 5.98 X 10 4 M for model V07.cs03. 



Table 2. Representative quantities of the considered models after quasi-stationary state has been reached. The columns report the model name, the average 
radiation temperature, the average effective adiabatic index, the accretion rate, the luminosity and the radiative efficiency, all of them computed after a quasi- 
stationarity state had been reached. See text for definition of these quantities. 



Name 


<T rad )[K] 


(7eff) 


M/A/ Edd 




'Ibu 


V07.CS03 


5.6 x 10 5 


1.333 


132 


0.939 


6.9 x 10" 3 


sp.V07.cs03 


5.6 x 10 5 


1.334 


135 


0.943 


6.8 x 10~ 3 


sp.V07.cs05 


4.3 x 10 5 


1.333 


62 


0.484 


9.0 x 10" 3 



ton fluid. 

Additional understanding of the thermodynamics of the models is 
achieved if we look at the time evolution of the averaged fluid and 
averaged radiation temperatures, which are plotted in Fig. [7] There 
are three points worth noting. First of all, for each model, the two 
temperatures T ra d and Tfluid differ by many orders of magnitude, 
suggesting that, at least globally, there is a strong deviation from 
thermal equilibrium within the fluid. Secondly, the fluid tempera- 
tures of the models sp.V07.cs03 and V07.cs03 are also signifi- 



cantly different, in spite of the dynamics being very similar (this is 
discussed in the next §). Finally, T ra d shows a smooth evolution, 
whereas T nu id exhibits a strong dip, reaches a minimum, and heats 
up again afterwards. When the large size andhot(T fluid > 10 10 K) 
shock cone reverses, the density downstream of the accretor be- 
comes small, yet the pressure remains high. A smaller size high 
temperature shock cone forms in the downstream region, as visible 
in the right panels of Fig. [6] with Tn u id-average being dominated 
by the high values within the shock cone. Thus, the fluid behaves 
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Figure 10. RADIATIVE EFFICIENCY - Comparison of the radiative effi- 
ciency f? BH£ as a function of time. 

like an effective photon fluid of temperature T ra( j<10 6 K, but in the 
shock cone no thermalization is possible and the fluid temperature 
vastly exceeds the radiation temperature. 

In order to corroborate this description, we measure the entropy 
generation rate V^S" as an effective tracer of dissipative processes 
(see AppendixfA}. The entropy generation rate ( Fig.[8]for the two 
models V07.cs03 and sp.V07.cs03) is maximum in the shocked 
region downstream of the accretor, where the fluid is very far from 
thermalization, even though the dynamics is otherwise stationary. 
Only in the upstream regions where V^S" is very small (cf. the 
white region of Fig. [8j, are the two temperatures T ra d and Ta u id 
similar. 



To further illustrate the effects of radiation induced dynamics, we 
measure two crucial parameters of accretion flows, namely the lu- 
minosity and the accretion rate, both of them reported in Eddington 
units. The computation is performed by directly integrating the es- 
caping radiation fluxes f r (Eq, IA7t . and the infalling mass fluxes 
at the innermost grid-point, respectively. We plot the time evolu- 
tion of the luminosity and of the accretion rate on the left and 
on the right panels of Fig. [9] for all three models. The important 
aspects of this figure are that (i) there is a transient peak in the 
luminosity evolution, corresponding to the point in the dynamics 
where the shock cone is momentarily dissipated away; (ii) the final 
luminosity is sub-Eddington for all cases; (iii) the luminosities of 
the models sp.V07.cs03 and V07.cs03 converge towards the same 
value; and (iv) the higher sound speed (correspondingly the lower 
asymptotic Mach number A4oo) of sp.V07.cs05 leads to smaller 
luminosity. On the other hand, the corresponding accretion rates are 
substantially super-Eddington, with final values of M / AfEdd in the 
range [62, 135], and confirming the advection dominated nature of 
BHL accretion flows. The relaxed luminosity efficiency rj BHC of 
the models together with all radiation quantities are listed in Tab.[2] 

The role of fluid temperature It is interesting to note that the 
strongly perturbed model sp.V07.cs03 converges towards a final 
state that is very similar to that of its unperturbed counterpart 
V07.cs03. This is observed in the accretion rate, M, plotted on the 
right panel of Fig. [9] in the luminosity, in the radiative efficiency, 77 
(cf. Fig.llOt. in the optical depth, r (cf. Fig. [6}, and in the radiation 
temperature, T ra d (cf. Fig. [Tj of these two models. 

This effect remained obscured in paperl, due to the fact that 



the previous criterion for the luminosity extraction was spuriously 
affected by boundary effectf^l 

While the radiation quantities converge for the two models 
sp.V07.cs03 and V07.cs03, the quantities more directly related to 
the fluid properties do not. For example, the fluid temperature, the 
Mach number and the entropy generation are neither qualitatively 
and certainly not quantitatively the same. In addition, the radia- 
tion dominated regime is reached at earlier times for V07.cs03 as 
it has higher Tfl u id and thus higher thermal conductivity (cf. Ap- 
pendix [A). 

From the consideration above, it stands to reason that the ra- 
diation temperature and the matter density (conversely, the optical 
depth) are the quantities affecting the dynamics most. This means 
that bremsstrahlung cannot not be a dominant process, since it is 
a temperature dependent radiation interaction. This is confirmed 
by the fact that, when looking at the respective opacities, Thom- 
son scattering dominates over bremsstrahlung by several orders of 
magnitude. We also note that in some portions of the grid, the dis- 
crepancy between Tfl u id and T ra d is very large, implying that the 
assumption of LTE is not valid there (cf. the red region of Fig. [8}. 
This is consistent with the fact that full thermalization in general 
is very hard to accomplish in dynamical environments of moderate 
density. 

Further comments We had already pointed out in paperl that 
models with initial high Mach number, Moo, are characterized 
by a luminosity that is dominated by the emission at the shock 
front, rather than by accretion-powered luminosity. This is also 
confirmed by the relative comparison between sp.V07.cs03 and 
sp. V07.cs05, the former having a larger Mach number and a higher 
luminosity (left panel of Fig.[9jl. 

Finally, we would like to comment about what has been 
dubbed the "'flip-flop'" instability in BHL accret ion flows, and 
whos e physical nature is still a matter of debate I Foglizz o et al.l 
2005). While we do not see this instability in our models (nei- 
ther in paperl nor in the present), we have observed that during 
the '"shock-reversal"' strong, although transient, oscillations in the 
shock cone can appear. However, we suspect that this effect can 
be partly attributed to the numerics, since the use of the IMEX 
in combination with a higher order Runge Kutta (order 3 instead 
of 2), alters the behaviour of this oscillation slightly. An extended 
analysis through three dimensional simulations would be needed to 
establish the potential relation of this oscillatory behaviour with the 
eventual development of the flip-flop instability. 



6 CONCLUDING REMARKS 

In this paper, we have revisited the optically thick, thermal radia- 
tion transfer in GR. First, we addressed the numerical problem of 
stiff source terms; proposed a numerical treatment, implemented 
and verified it. As we chose an IMEX Runge-Kutta scheme, we 
needed to isolate the principal stiff parameters, which were found to 
be the (density-weighted) opacities. After applying the new IMEX 
method to the one-dimensional problem of spherical a ccretion, we 
compa red our results with those obtained earlier by No bili et all 
( 1991) and found good agreement. In this spherical, stationary sce- 
nario the current formulation of the GR-RHD equations is fully ap- 
plicable as long as the solution remains optically thick. We remark 

14 See discussion in Appendix lAl 
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Figure 8. 2D Entropy generation rates for V07.cs03 and sp.V07.cs03 - 
Distribution map of V^S" for model V07.cs03 (left panel) at time t = 5.98 X 10 4 M and for model sp.V07.cs03 (right panel) at time t = 7.71 X 10 4 M. 
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Figure 9. Time evolution of perturbed BHL models sp.V07.cs05, sp.V07.cs03 and V07.cs03- 
(Left panel) luminosity L extracted at constant optical depth r > 10, (Right panel) accretion rate M as a function of time in Eddington units. 



that there is not a unique stiffness threshold, valid for any physical 
scenario, at which the purely explicit scheme fails and the IMEX 
becomes necessary. In the case of a purely explicit RK scheme, 
when the source terms become stiff, it is possible to a certain extent 
to lower the Ccfl factor and obtain a stable evolution. However, 
the stiffness parameter can become very large, so the time-step very 
tiny. That is of course inefficient, and resorting to a stiff solver is 
the only way out. In general, if a problem can be solved with a 
purely explicit RK scheme, this is to be preferred as it is CPU- 
faster. However, we believe that most non-trivial radiation applica- 
tions will exhibit stiffness and lead to code crashes with standard 
explicit RK schemes. 



We then revisited the Bondi Hoyle Lyttleton accretion in 2D 
for an astrophysical, dynamical problem. Here, we could show that: 

• The IMEX scheme allows us to evolve models with realistic 
choice of parameters, of the order of T ~ 10 6 K; 

• the dynamics of the flow are significantly affected by the ra- 
diation pressure, yielding super-Eddington accretion rates in the 
range M ~ [62, 135]MEdd and Eddington limited luminosities; 

• the fluid and the radiation depart strongly from thermal equi- 
librium in shocked regions, particularly in the shock cone down- 
stream of the accretor. 

Our analysis has substantially benefited from the ability of our 
scheme to treat stiff source terms. However, we should also state 
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a few words of caution as to the current shortcomings and neces- 
sary future improvements of our scheme: 

• The optically thin regime cannot be treated yet, and further 
steps are required to incorporate the variable-Eddington factor ap- 
proach. 

• Temperatures of order T < 10 K, as they appear in small re- 
gions of the domain, require the inclusion of bound-free opacities, 
which are currently neglected. 

• The only dissipative mechanism is currently thermal conduc- 
tivity. Other types of viscosity such as an effective viscosity related 
to magnetic turbulence would be beneficial. Coupling the current 
equations to MHD represents another direction of future research. 

• Since we currently cannot extract the luminosity in regions 
where the optical depth is low, we must trace a geometrical surface 
of constant r ^ 1. However, it remains an uncertainty as to where 
such a surface should be placed, and the computed luminosities are 
therefore affected by at least one order of magnitude uncertainty. 

Even in the presence of these limitations, our analysis may become 
relevant for the study of merging supermassive black-hole bina- 
ries, which have been attracting a lot of interest for the possible 
joint measure of electromagnetic and gravitational wave signal (in 
the context of multi-messenger a stronomy). Neglec ting the back- 
reaction of radiation onto matter, iFarrisetalJfcOlOh already con- 
sidered the BHL solution in a binary system, finding that luminosi- 
ties as high as 10 43 ergs _1 can be obtained in a hot gas cloud of 
temperatures T ~ 10 6 K. Such estimates are compatible with our 
calculations, but a dedicated work will be presented in the future. 
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where a M is the four-acceleration of the fluid, A is the thermal con- 
ductivity and h^v — g M „ + u^Uv is the projector operator in the 
space orthogonal to the four- velocity tt M . Under the assumption 
that the fluid has vanishing shear and vanishing bulk viscosity, the 
entropy-generation rate that follows from JAU and JA2t is given 
by 



AT 



(A3) 



We recall that the thermal conductivity is related to the opacity. 
For instance, the thermal conductivity computed using the ordi- 
nary diffusion app roximation of st ellar interiors is given by A = 
(4/3)a rad cT 3 /x s dSchwartJl967h . Under the assumption that the 
matter plus radiation fluid behaves as a single fluid with effective 
pressure and energy density given by P e g = P+Vi, e e g = e + E r , 
the four acceleration a M can be computed from the Euler equations 
as 



(A4) 



e cB + P eB 

When quasi stationary configurations are reached, the terms con- 
taining time derivatives can be neglected with respect to those con- 
taining spatial derivatives, and after replacing q M into Eq. iA3\ we 
obtain 



J"! 



( 9 rr + rV) xa-rf + 



(g™ + r 2 (t/) 2 )(9 T) 2 + 2T 2 v r v'*'drTd4,T - 



2T 



e c ff + PeB 



+ T 2 (v 4 



(g rr + r 2 (v r ) 2 )d r P aB drT + 

2 )d^Td^P cS + 
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T 



e c ff + Pes 
(g^ + r 2 (v^) 2 )(d 4 ,P aH ) 2 + 

2T 2 v r v 4 'drP a sd 4 ,P e s 



+ Y 2 {v r ) 2 ){drP, a f + 



(A5) 



The conversion of V M 5 M from geometrized units to cgs units is 
given by 

3 



[V^5 M ] cgs = 1.0353xlO J1 Gc[-jZ] [V M 5 



KL 



(A6) 



In the code we generally compute the luminosity as the surface 
integral over outgoing radiation fluxes f£ as 



APPENDIX A: ENTROPY GENERATION RATE AND 
LUMINOSITY COMPUTATION 

In the framework of Eckart's formulation o f Relativistic Standard 
Irreversible Thermodynamics teckartlll940h . the entropy current is 
given by 



(Al) 



where g M is the heat flux, s is the entropy per unit mass, and T is 
the temperature of the fluid. The heat flux is given by the relativistic 
form of Fourier law, namely dlsraellll976h 



g M = -AT(XV„ In T + a» 



(A2) 



L = 2 £ [V7 (/;)„ A0„ 



(A7) 



where A<fi n is the angular size of a grid cell and the integral is taken 
at the radial position of the last optically-thick surfaca 15 L i.e. where 
t — t.. In paperl we computed the luminosity by imposing the 
criterion r. ^ 1. However, these small values of the optical depth 
often correspond to an integration surface close to the boundary of 
the numerical domain, where spurious boundary effects may alter 
the results. Hence, in this paper we have adopted a different cri- 
terion by choosing r. ^ 10, which guarantees that the integration 

15 The factor 2 in jA7t accounts for both the contributions above and below 
the equatorial plane. 
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Figure Al. Comparison of luminosity extraction of perturbed 
BHL - p.V18.cs07 and p.V10.cs07. Extracting the luminosity at r ^ 10 
leads to different light curves, these curves are labelled "new". In paperl we 
had used the criterion r > 1. 
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Figure A2. Uncertainties and comparison of luminosity ex- 
traction - V09.cs07: Extracting the luminosity at r ^ 10 leads to differ- 
ent light curves with much smaller uncertainties, these curves are labelled 
"new". In paperl we had used the criterion t > 1. 



surface is not placed at the outermost grid cells. For clarification we 
have repeated the luminosity extraction for two models considered 
in paperl, p.V18.cs07 and p.V10.cs07, and show the light curves, 
computed with the two different criteria, in Fig. lAll 

We can assign error bars to our extraction method by taking 
the standard deviation of the mean. For model p.V09.cs07, the 
comparison is shown, including the errorbars, in Fig. lA2l We stress 
that the size of such error bars reflects the uncertainty in choos- 
ing the position of the last optically thick surface across which the 
emitted luminosity is computed. It should be noted, moreover, that 
both our estimates agree within these uncertainties, but the choice 
r ^ 10 produces much smaller error bars than r ^ 1 and should 
therefore be preferred. 



APPENDIX B: IMPLEMENTATION OF THE IMEX 
SCHEME 

A tableau notation is usually adopted to express in a compact form 
the coefficients of the matrices a,ij, Uij and of the corresponding 
vectors u)i, Cbi as 



(Bl) 



where the index T denotes transposition^ 

The explicit tableau of the SSP3(4, 3, 3) is 
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; 0.24169426078821 
0.12915286960590 



q 2 = 0.06042356519705 , 



The coefficients of the radiation matrix A(Y) of Eq. d37b are 
given by 
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16 Note that the coefficients and which are defined as c; = 
2~Zj = i o-ij and iii = X/j=i <Hj ^ n °t us ^d in the practical implemen- 
tation of the scheme. 
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where, just for convenience, we have specified the spatial coordi- 
nates to (x, y, z). 
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